Bond disorder induced criticality of the three-color Ashkin- Teller model 
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An intriguing result of statistical mechanics is that a first-order phase transition can be rounded by disorder 
coupled to energy-like variables. In fact, even more intriguing is that the rounding may manifest itself as a 
critical point, quantum or classical. In general, it is not known, however, what universality classes, if any, 
such criticalities belong to. In order to shed light on this question we examine in detail the disordered three- 
color Ashkin-Teller model by Monte Carlo methods. Extensive analyses indicate that the critical exponents 
define a new universality class. We show that the rounding of the first-order transition of the pure model 
due to the impurities is manifested as criticality. However, the magnetization critical exponent, /3, and the 
correlation length critical exponent, v, are found to vary with disorder and the four-spin coupling strength, and 
we conclusively rule out that the model belongs to the universality class of the two-dimensional Ising model. 

PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q, 64.60.-i 



Rounding of first-order transitions due to disorder is an im- 
portant problem because such transitions are ubiquitous in 
both classical and quantum systems. In contrast to the ef- 
fect of disorder on continuous transitions, much less is known 
about its effect on first-order transitions. Imry and Wortis flj, 
Hui and Berker |2j, and Aizenman and Wehr [3] have made 
important contributions to this subject. 

More recently, the quantum A-color Ashkin-Teller (AT) 
model in one dimension was studied using both weak and 
strong disorder analyses |4| where evidence was presented 
in favor of rounding of this transition into a quantum criti- 
cal point for N > 3. An important question that is unan- 
swered is the universality class of this transition. It was ar- 
gued, however, that the physical picture is the same as the 
random transverse-field Ising model for a limited parameter 
regime. Even more recently, general and rigorous analyses 
have shown that disorder can indeed round a first-order quan- 
tum phase transition |5 1. However, it is still unknown what the 
possible universality classes are, or if rounding is not simply 
a smearing of the first-order transition instead of an emergent 
criticality. 

The phase transition in the pure iV-color AT model has been 
studied for decades @. It is known that for A > 3, the pure 
model undergoes a first-order phase transition in two dimen- 
sions {2D). However, the corresponding disordered case is 
still far from being understood. Here, we address the criti- 
cality of the bond disordered classical three-color AT model 
in 2D, a necessary prerequisite for understanding the corre- 
sponding quantum phase transition at zero temperature. Uti- 
lizing the Monte Carlo method, we compute the critical ex- 
ponents for the magnetization and the localization length, and 
show that the phase transition does not belong to the Ising uni- 
versality class, as previously suggested Q. We also find that 
the critical exponents vary with disorder and the four spin cou- 
pling strength. Previously a perturbative two-loop renormal- 
ization group calculation had hinted at a new strong-disorder 
fixed point [8]. That the Ising universality class does not hold 



appears to be similar to the g-state Potts model |9| in the pres- 
ence of disorder, where the exponent also varies with q. 
The Hamiltonian of the A-color AT model in 2D is 
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where the classical spin variables sf = ±1 reside on a square 
lattice, and the N spins on a single lattice site are labeled by 
a, (3 = 1,2, . . . , N, The sum is only over nearest neighbor 
pairs. In our study, N = 3 and we only consider the case 
for g > 0. The ferromagnetic coupling constant Jjj > 
is a random variable following a translation-invariant binary 
probability distribution of the form: 
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The magnetization of the system is 



1 

N 



I N 

E 



(3) 



where the symmetry between spins of different colors is uti- 
lized to increase accuracy. The angular brackets (• ■ ■ ) denote 
the usual thermal Monte Carlo average, whereas the square 
brackets [• ■ ■ ] denote the quenched average over configura- 
tions with different {Jij}- 

Close to a continuous transition, the magnetization m 
scales as 
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where rh(-) is a universal function, and v is the critical expo- 
nent for the correlation length. The variable x = ( J — J c ) / J c 
is the reduced coupling constant in temperature units. 

There are two other quantities of interest that are useful in 
determining the nature of the phase transitions: the energy 
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cumulant V c and magnetic cumulant V m ifTUl [TT l given by 
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respectively, where e = T-L/(NL 2 ), H is the Hamiltonian 
given by Eq. ([TJ, and L is the linear dimension of the square 
lattice. For a continuous transition, the cumulant V m has the 
scaling form 



V rn = V m (XL 1 '"} 
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Using V e we can determine the order of the phase transition. 
In a continuous phase transition, V e — > 2/3 V J as L — > oo, 
whereas in a first-order phase transition, V e behaves like the 
continuous case except at J = J c , where approaches a 
non-universal constant. 

To apply the cluster Monte Carlo algorithm to our system, 
we fix a color A and rewrite the Hamiltonian in Eq. ([T| to 
obtain 



nonzero finite value. This validates that the transitions are 
all first-order for g < 0.18. Moreover, expecting the transi- 
tion point to be close to the minimum of V £ for the largest 
examined size, we deduce the phase diagram of the pure AT 
Hamiltonian as shown in Fig.[T] 
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where iJj represents the terms which do not include spins of 
color A, and Efj is the bond energy for the bond between 
sites i and j in color A. This expression can be viewed as 
a random Ising model Hamiltonian with equivalent nearest- 
neighbor coupling constant Jj j + g Y2 a £A s i s< j ^ oni y c °l° r 
A is considered. Therefore, any cluster algorithms for the ran- 
dom Ising model can be adopted here. For the following cal- 
culation, the algorithm proposed by Niedermayer [ 12 1 is used. 
Note that the bond energy is bounded from above by 



E max = J + — + (N — l)g. 
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The equilibration time for all observables are found by log- 
arithmic binning. We perform one million Monte Carlo steps 
for equilibration. Furthermore, we perform 10, 000 thermal 
averages over each disordered configuration. The number of 
disorder averages used for each observable varies and is given 
in the caption of the plots. 

Because all temperatures are simulated with the same dis- 
order realization, the measured data are correlated. Therefore, 
in order to obtain the error for the cumulants, we have applied 
the Jacknife procedure to correct for bias JT3). 

First, we must identify the parameter regime in which the 
pure system undergoes a first-order transition by looking at 
the phase space of the pure AT Hamiltonian. For fixed val- 
ues of g, we calculate V t for system sizes L = 24, 32, . . . , 64. 
The first-order transition is confirmed by inspecting how the 
depth of the energy cumulant, 2/3 — V*[J C (L), L], changes 
as a function of some inverse power of the system size. As 
expected, the extrapolated lines intersect the ordinate at some 



FIG. 1. (Color online) Phase diagram of g versus J with A = for 
TV = 3. For any values of g < 0.18, the system undergoes a first- 
order phase transition. The error bars are smaller than the size of the 
symbols used in the plot. 

Now we study the disordered system for the parameter set 
(g, A) = (0.1, 0.2) and use V e , shown in Fig. [2] to find the or- 
der of the transition. Plotting the depth of V t against 1 /L 3 / 4 
(Fig. [3]), we observe that V* — > 2/3 as L — s- 00, which indi- 
cates that the transition is continuous. 
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FIG. 2. (Color online) The energy cumulant V e versus J with 
(g, A) = (0.1, 0.2) for TV = 3. The measured values are averaged 
over 15, 000 configurations. V e approaches the value 2/3 for all J in 
the continuous phase transition case and eventually becomes 2/3 for 
infinite system size. The error bars are hardly visible because they 
are smaller than the symbols. 

We use the scaling properties of V m in Eq. |6]), and the log- 
arithmic derivative of the magnetization squared to find an ac- 
curate estimate of the transition point J c . Figures |4]and [5] 
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FIG. 3. (Color online) (2/3 - V*(J C (L),L)) versus 1/L S/4 with 
(g, A) = (0.1,0.2) for N = 3. The line goes through the origin 
indicating a second-order phase transition. The error bars are smaller 
than the size of the symbols used in the plot. 
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FIG. 5. (Color online) Logarithmic derivative of magnetization 
dlog(m)/dJ versus J with (g, A) = (0.1,0.2) for N = 3. The 
measured values are averaged over 20, 000 configurations. The in- 
tersection point is at J c = 0.2328 ± 0.0003. The error bars are 
hardly visible because they are smaller than the symbols. 
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show that the two quantities are independent of the system 
size at a single critical point, J c = 0.2331 ± 0.0006. 
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FIG. 4. (Color online) Magnetic cumulant V m versus J with 
(g, A) = (0.1,0.2) for N — 3. The measured values are av- 
eraged over 49, 000 configurations. The intersection point is at 
J c = 0.2334 ± 0.0003. The error bars are hardly visible because 
they are smaller than the symbols. 

From the scaling behavior of V m in Eq. (l6t, we deduce the 
critical exponent v. The best data collapse is obtained with 
J c = 0.2335 ± 0.0001 and v = 0.76 ± 0.05, see Fig. [6] 
Furthermore, the critical exponents /3 and v can be extracted 
from Eq. (RJ. As shown in Fig. [7] the magnetization scales 
well with v = 0.70 ± 0.02 and f3 = 0.055 ± 0.005. 

To address the universality of the exponents, it is important 
to study the behavior of the system for more than just one 
particular set of parameters (g, A). In this regard, we repeat 
our calculation for the parameter sets (g, A) = (0.05, 0.20) 
and (0.10,0.10). Table [I] summarizes the values of v and /3 
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FIG. 6. (Color online) Finite size scaling plot of the magnetic cu- 
mulant with (g, A) = (0.1,0.2) for N = 3, J c = 0.2335, and 
v — 0.76. The upper branch corresponds to the ordered phase 
whereas the lower branch corresponds to the disordered phase. The 
error bars are hardly visible because they are smaller than the sym- 
bols. 



for each parameter set. 

Based on our results, the disordered three-color AT model 
does not belong to the Ising universality class. One is tempted 
to conclude that reducing the coupling constant g has no effect 
on the exponent v (at least within the error bars) but results in 
an increase in the exponent /3, On the other hand, reducing 
the disorder A has a lesser impact on the exponent j3, but it 
results in smaller v. 

The results, when taken at face value, seem to violate the 
bound v > 2/D lfl4ll . The value of v is extremely sensitive 
to J c : it is remarkable that a variation of J c from 0.2334 to 
0.2336 can change v from 0.71 to 0.81. In contrast, the same 
variation of J c results in a much weaker variation of the ex- 
ponent j3 from Eq. (HI). The spread of (3 is in the range of 0.05 
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FIG. 7. (Color online) Finite size scaling plot of the magnetiza- 
tion with (g,A) = (0.1,0.2) for TV = 3, J c = 0.2334, and 
/3/v « 0.07. The magnetization values are averaged over 10,000 
configurations. The upper branch corresponds to the ordered phase 
whereas the lower branch corresponds to the disordered phase. The 
error bars are hardly visible because they are smaller than the sym- 
bols. 





(0.10,0.20) 


(0.05,0.20) 


(0.10,0.10) 


v m 


v = 0.76 ±0.05 


i/ = 0.71 ±0.01 


v = 0.59 ±0.04 


m 


v = 0.70 ±0.02 
,5 = 0.055 ±0.005 


v = 0.79 ±0.02 
f3 = 0.080 ±0.005 


v = 0.53 ±0.01 
p = 0.060 ± 0.007 



TABLE I. The values of critical exponents extracted by finite-size 
scaling of magnetic cumulant and magnetization for different values 
of coupling constants and disorder. 




FIG. 8. (Color online) Finite size scaling plot of the magnetic cu- 
mulant with (g, A) = (0.1,0.1) for N = 3, J c = 0.2436, and 
v = 0.56. The upper branch corresponds to the ordered phase 
whereas the lower branch corresponds to the disordered phase. The 
error bars are hardly visible because they are smaller than the sym- 
bols. 



to 0.06 which is considerably smaller than the Ising exponent 
1/8. 



Another striking fact is that the critical exponents (3 and v 
vary as the coupling parameter g or the disordered strength 
A is altered. Similar behavior has been previously reported in 
studies of the critical behavior of the ^-states Potts model [15], 
where the magnetic exponent varies continuously with q > 4. 
Interestingly, in the very same studies the bound v > 2/D 
was violated, as in our present analysis. 

It is worth mentioning that we have also examined the 
finite-size scaling of the correlation length £ extracted from 
the second moment correlation function |[T6l 
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The correlation length scales as £/L = £ (xL 1 /") for a con- 
tinuous transition. Although preliminary results for v are sim- 
ilar, the larger system sizes are much harder to achieve with 
our limited computation time, since the sum in Eq. ( fT0| ) is to 
be performed over all lattice points. The multifractality of the 
correlation function 1 9 1 may be another difficulty in extracting 
reliable results. 

In summary, although precise values of the exponents are 
difficult to establish, the emergence of criticality induced by 
disorder in the pure AT model for N = 3 (and presumably 
for N > 3 as well) is manifest and the presence of a sin- 
gle universality class (such as Ising) appears to be excluded. 
Further analytical and numerical studies of the corresponding 
quantum version [4] remain to be conducted. 
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